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Abstract 

We propose and implement a new method, called the Ice Sheet Cou¬ 
pled Approximation Levels (ISCAL) method, for simulation of ice sheet 
flow in large domains under long time-intervals. The method couples 
the exact, full Stokes (FS) equations with the Shallow Ice Approximation 
(SIA). The part of the domain where SI A is applied is determined auto¬ 
matically and dynamically based on estimates of the modeling error. For 
a three dimensional model problem where the number of degrees of free¬ 
dom is comparable to a real world application, ISCAL performs almost 
an order of magnitude faster with a low reduction in accuracy compared 
to a monolithic FS. Furthermore, ISCAL is shown to be able to detect 
rapid dynamic changes in the flow. Three different error estimations are 
applied and compared. Finally, ISCAL is applied to the Greenland Ice 
Sheet, proving ISCAL to be a potential valuable tool for the ice sheet 
modeling community. 


1 Introduction 

With the increasing awareness of climate change and the associated sea level 
rise, research on ice sheets has intensified [3, 18, 31, 37, 40]. Accurate predic- 
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tions of the response of the Greenland Ice Sheet and the Antarctic Ice Sheet 
to global warming call for accurate numerical models, with an adequate repre¬ 
sentation of ice dynamics. The most accurate models of ice dynamics available 
to date, are so-called Full Stokes models (hereafter referred to as FS). In these 
models the ice dynamics are represented by the Stokes equations, treating ice 
as a slowly creeping, non-Newtonian, viscous fluid [17]. These equations are 
computationally demanding to solve since they are highly non-linear, with a 
viscosity which depends on both strain-rate and temperature. 

Computer performance restricts ice sheet simulations with FS to rather short 
timescales on the order of a few hundred years even on the largest parallel com¬ 
puters, see [35, 16, 25, 30, 29, 42]. Yet, the overall very slow response of ice 
sheets to external forcings (typical response times are thousands of years) calls 
for investigations of ice dynamics on much longer timescales. Longer simula¬ 
tions allows us to determine whether observed short-term fluctuations represent 
natural variability or are part of a long-term trend. Furthermore, we can vali¬ 
date ice sheet models by comparison to geological data, which provide a record 
of the long-term evolution of ice sheets that took place before the beginning 
of the modern, mostly satellite-based, observational period. Such simulations 
of long-term ice sheet behavior are known as paleo-simulations, which typically 
also involve considerably larger spatial domains than those currently covered. 
These simulations are not yet feasible with FS. 

Numerical simulations of ice sheet behavior during glacial cycles at hemi¬ 
spheric scales are possible [38, 8, 39, 12] if an approximation to the Stokes 
equations is used, known as the Shallow Ice Approximation (hereafter referred 
to as SIA) [5, 13, 19]. The SIA is based on the assumption that the ice body 
is shallow, and that vertical shearing motion dominates over horizontal stretch¬ 
ing and shearing. This implies substantial simplifications to the FS equations, 
thereby rendering SIA models computationally very cheap. 

The overall, inland, ice sheet behavior is often modeled sufficiently well with 
SIA, but is poorly represented in regions which are now regarded as critical for 
the mass budget of ice sheets and for future sea level rise. Such regions are 
the highly dynamic ice sheet margin, where ice streams (faster flowing regions 
within the ice) terminate at the ocean. In the ocean, the ice either breaks into 
icebergs, or forms large floating ice platforms which remain connected to the ice 
sheet, called ice shelves. At ice shelves, and at domes or over rapidly undulating 
basal terrain, an accurate representation of ice dynamics is also not possible in 
SIA [5, 34, 20, I, 2]. As a consequence of approximations being applicable in 
some regions of an ice sheet, and in some other not, the idea of coupling the FS 
with approximations emerged in the ice sheet modeling community, presented 
first by Seroussi et al [36] in 2012. In [36], a coupling between the so-called 
Blatter-Pattyn approximation [6, 28], the Shelfy Stream Approximation [27], 
and FS was demonstrated for an idealized floating ice shelf. 

Theoretical analyses and model case studies have been devoted to understand 
precisely under what conditions and in which parts of an ice sheet SIA can be 
used and is expected to give accurate results [19, 34, 1, 26, 14]. Yet, it is not easy 
to transfer such insights to simulations of real ice sheets, featuring complicated 
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geometry over complex topography, spatio-temporally variable basal conditions 
and climate forcings. It is practically impossible to determine a priori when, 
and in which parts of a sheet, the SIA is applicable. 

Here, we propose a new method, called the Ice Sheet Coupled Approxima¬ 
tion Level (hereafter referred to as ISCAL). This method automatically and 
dynamically decides from within the simulation where in an ice sheet the SIA 
is a valid approximation, and applies the FS elsewhere. The ISCAL method is 
computationally less expensive than FS, and more accurate than the SIA. The 
ISCAL method goes beyond the method in [36], providing an automatic switch 
between FS and SIA based on estimated modeling errors, and coupling the FS 
and the SIA in a dynamic manner rather than statically. ISCAL has the poten¬ 
tial to simulate ice sheet dynamics on larger spatial and temporal domains that 
were previously considered inaccessible with FS alone. This opens the way to 
more accurate paleo-simulations, model validation by comparison to geological 
archives, and improved assessment of observed ice dynamics. ISCAL can also 
speed up simulations on shorter time intervals for parameter sensitivity studies 
or spin-up simulations to determine consistent initial conditions for FS where 
the accuracy requirements are lower. 

The aim of this paper is to present the ISCAL algorithm, which is imple¬ 
mented in the finite element software Elmer/Ice [15], and to evaluate its proper¬ 
ties. In Section 2, the theoretical framework of FS and SIA is briefly reviewed. 
In Section 3, the dynamic coupling of FS and SIA, and the procedureof how to 
estimate where the SIA is valid to a given accuracy, is described. In Section 4, 
we perform numerical experiments to demonstrate the accuracy, efficiency and 
flexibility of ISCAL using a transient conceptual model problem in three dimen¬ 
sions (3D). We also apply ISCAL to real data from the Greenland Ice Sheet to 
prove the maturity of the method and its potential value to the glaciological 
community. 


2 Theory of FS and SIA 

On time scales relevant for ice sheet dynamics, ice can be described as an 
isotropic, viscous, incompressible, non-Newtonian fluid, with a very low Reynolds 
number. The dynamics is governed by the non-linear Stokes equations, i.e. the 
FS equations. The SIA is derived from the FS equations by assuming that the 
aspect ratio, i.e. the ratio e = [iL]/[I/] of a typical thickness of the ice [H] to a 
typical lateral length scale [L] is small [5, 1,2]. By expanding the solution in the 
small parameter and solving for the lowest order terms, the SIA is obtained (see 
[5, 20]). In this section, the FS and SIA equations are described. We consider 
isothermal conditions for simplicity, although ISCAL is implemented for a fully 
thermo-dynamically coupled setting. A partial derivative of a variable u with 
respect to an independent variable x is written d^u. 

The general setting is as follows: We consider a grounded ice sheet in a 
Cartesian coordinate system, see Fig. I, where the ice sheet surface 2 ; = h{x, y, t) 
is a function of position {x, y) and time t, and the ice sheet base is z = b{x, y, t). 
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Figure 1: An ice-sheet in a Cartesian coordinate system, exaggerated in the 
z-direction. The ice surface position is described by h(x,y,t) and the bedrock 
underneath is given by b{x,y,t). The normal vector, n is pointing outwards 
from the ice. In two dimensions there is one tangential vector t, and in 3D we 
choose two that span a tangential plane. 


The ice sheet surface evolves as 

dth -\- Vx\z—h(^xb^ T '^y\z—hdyh Vz — ; (1) 

where Ug = as{x, y, t) is a prescribed source term denoting the net accumulation 
normal to the surface of the ice sheet, and depends on the climate conditions. 
The velocity field v = {vx,Vy,Vz)'^ is obtained as the solution to the FS or SI A 
equations. 

The bedrock underneath the ice is here assumed to be rigid, so that b{x, y, t) = 
b(x,y). For simulations covering long time spans, an evolution equation for b 
similar to (1) is solved to account for isostatic adjustment and bottom melt¬ 
ing/refreezing. At the boundaries (at i.e. the ice surface and at the ice sheet 
base), a normal unit vector n points away from the ice body, and two vectors 
ti and t 2 span a tangential plane. 

2.1 The full Stokes formulation 

The FS equations that are solved for the velocity v = {vx,Vy,Vz)’^ and the 
pressure p are 


-Vp -I- V • (r 7 (Vv -f (Vv)”^)) -I- pg = 0, 

V • V = 0, 


(2a) 

(2b) 
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where p is the density, and pg is the force of gravity. The viscosity, ry, depends 
on the velocity as 




(3) 


The Glen exponent n is usually chosen as n = 3. The rate factor A = A{T') 
is coupling the pressure melting point corrected temperature field, T', to the 
velocity field through an exponential function [17]. Due to convection in the 
temperature equation, there is thus a two-way coupling between temperature 
and velocity. As we only consider isothermal conditions, we will not describe 
this here, and the value of A is simply a constant given in Section 4.1. The 
effective strain rate, d, is given by 


d = 


\/(ltr(5(Vv-H(Vv)^)^)) = ^ A {dyVyf A 

+ I {dyv^ A d^Vyf + i {d^v^ + A i {dvy + dyV^'^^ ^ ■ 


(4) 


Since typically n > 1, there is a singularity obtained in (3) at c? = 0, which 
can cause numerical difficulties. To regularize the problem, the shear rate is 
bounded below by a small number, do [32]. 

Due to the velocity dependent viscosity in (3), the problem (2) is highly 
non-linear. Becasue of this non-linearity, and due to the fact that the equation 
system in (2) is a saddle-point problem, the FS equations are computationally 
expensive to solve. 


2.2 Boundary conditions 

At the ice surface, the atmospheric stresses are generally considered to be neg¬ 
ligible implying 

cr • n = 0 (5) 

The Cauchy stress cr is related to the deviatoric stress tensor t, velocity and 
pressure by 

cr =-pi -I- X =-pi -I- ? 7 (Vv -I- (Vv)^), (6) 

i.e., at the surface we have a homogeneous Neumann condition for the velocity 
field and a homogeneous Dirichlet condition for the pressure. At the ice-bedrock 
interface, it is less obvious how to choose boundary conditions. The conditions 
there depend on geological, hydrological and thermal conditions which are often 
inaccessible to observations. We employ two different boundary conditions for 
the velocity v at the base (indicated by subscript b): a no slip condition 

v|6 = 0, (7) 

which is accurate if the ice base is frozen, or a linear sliding law, relating the 
basal velocities with basal shear stress as 

(v • = -Th//3 = -(t; • cr • n)|t,//3, i = 1,2 

(v • n)\b = 0. 


5 


(8a) 

(8b) 



The vectors ti and t 2 span a tangential plane, and the normal vector n is 
pointing downwards (Fig. 1). We assume that no melting or refreezing occurs 
at the bed, hence (8b) holds. The friction between the ice and the underlying 
bedrock is described by the sliding coefficient /3, which will be small in high 
sliding areas such as under ice streams/outlet and glaciers. 


2.3 The Shallow Ice Approximation 

Due to the approximations made in deriving SIA, the SIA velocities and pressure 
are very simple to compute and are directly given by: 


Vx 


V 


y 


= t'M-2(P5)”a,h||v/i|r-i 

= Vb,y-2{pgrdyh\\^hr-^ 


r A{T'){h-z'Tdz\ 

Jb 

f A{T'){h-z'Tdz\ 
Jb 


Vz = Vb,Z - / {dxVx + dyVy) Jz' , 

Jb 

P = P9{h-z), 


(9a) 

(9b) 

(9c) 

(9d) 


see [17]. The norms in (9) and in all following equations is the Euclidean vector 
norm. To arrive at the integral forms of (9), the boundary conditions were 
used, and except for a sliding law and a melting/refreezing model describing 
the basal velocity {vb,x,Vb,y,Vb,z), no further equations are needed. The sliding 
law is the linear sliding law in (8), but the Cauchy stress tensor at the base, cr, 
is approximated as 


/ -pg(h-b) 0 -pgd^h{h-b) \ 

(T = I 0 -pg{h-b) -pgdyh{h-b) J . (10) 

\-pgd^h{h-b) -pgdyh{h-b) -pg{h-b) / 

As in (8b), we assume there is no melting or refreezing at the ice base. 

We previously examined the validity and accuracy of the SIA in [1] and [2] 
by analysis and numerical simulations of 2D ice flow over a bumpy bed. As 
expected, the accuracy of the SIA decreases when the aspect ratio e increases. 
For aspect ratios above 10“^ the relative error is of 0(1). Furthermore, the SIA 
theory assumes that the surface slope is proportional to £, [2], that shearing 
dominates sliding, and that horizontal in relation to vertical velocities scale 
with the aspect ratio e [34]. The SIA is thus not expected to perform well in 
regions with large spatial variations in data (e.g. topography), at steep margins, 
in ice streams or ice shelves, and at domes. 

2.4 Time evolution with FS and SIA 

The FS and SIA equations (2) and (9) are time-independent. Time-dependency 
enters in the free surface problem in (I), and in the temperature equation (not 
considered here), which in turn determines the evolution of the viscosity. The 
surface evolution equation (I) exists in both a non-approximated version (1) 
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Figure 2: An extruded, 3D mesh on a circular ice sheet, with 17860 nodes, 
distributed in 20 vertical layers. The vertical axis in the graph is scaled by a 
factor 100 in the figure. A slice is cut out to show how the nodes are aligned in 
the vertical direction. This 2D footprint mesh is unstructured in the center and 
structured at the margins. 


and an approximated SIA-version. The computational gain in approximating 
them with SIA is far less than for the solution of the Stokes equations. In our 
simulations, we stick to the solver already implemented in Elmer/Ice [15]. The 
following discussion concerning the ISCAL method will focus on the computa¬ 
tion of the velocity and pressure fields. 

3 Numerical method 

The numerical solution of the ice sheet evolution consists of two parts. First, 
the FS and SIA equations (2) and (9) are solved for the velocity v with a fixed 
ice domain and a computational mesh covering the domain at time t. Then a 
new height at t -I- At is computed by integrating (1) in time using the velocity 
computed at t. The mesh is then adjusted in the vertical direction to account for 
the elevation changes and (2) and (9) are solved again at t -|- At. The backward 
Euler time integration is first order accurate with an error of 0{At). 

Since ISCAL is based on the FS and SIA equations, we first briefly describe 
how to solve the respective equations and then proceed to the combination in 
the ISCAL method in more detail. 

3.1 Solving the full Stokes equations 

Finite element discretizations of the FS equations, the ice surface evolution, 
and the temperature evolution are already implemented in Elmer/Ice [15]. The 
mesh consists of triangles or quadrilaterals in the x — y plane extruded in the 
z direction to the ice surface to form prisms, see Fig. 2. This type of mesh is 
favorable for ice sheet modeling. Linear elements are used, and the discretization 
of (2) with a given (e.g. from a previous iteration) viscosity, y, leads to a system 
of linear equations 

Au = b, (11) 
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where the components of the solution vector u are the nodal values of v and 
p. The right hand side, b, is an expression of the gravitational force for the 
momentum equations. The solution of (11) with a fixed A is computed by 
an iterative method. Because of the non-linear rheology of ice, the system 
matrix A depends on 77 which in turn depends on the velocity, v, through 
(3). Hence, A = A(r7(u)) = A(u). Either a Picard iteration or a Newton 
iteration is applied to solve the non-linear equation (11), generating a sequence 
of solutions Ufe, ( k = 1,2,...). Both the linear system iterations and the non¬ 
linear iterations are considered as converged when the difference (measured in 
II • II) between two consecutive and u^, is sufficiently small. The initial 
guess Uq in the non-linear iterations is the solution from the previous time step. 
This iteration in combination with large domains and long time intervals are the 
main reasons why ice sheet modeling is a computationally demanding problem. 

3.2 Solving the SIA equations 

The SIA solution is determined by computing the right hand side expressions 
in (9). As there are no partial differential equations to solve, the numerics is 
simple. The gradients in (9) are computed using the finite element framework 
of Elmer/Ice. As the meshes are constructed such that the nodes are aligned in 
the vertical direction, (Fig. 2), the integrals in (9) can be computed using the 
trapezoidal rule. 


3.3 Coupling FS and SIA 

Let us introduce some specific notation for describing the coupling. The SIA 
solution is denoted by ug/^, the FS solution by and we let the combined 
ISCAL solution, obtained by coupling FS and SIA, be written ujscal- The 
computational domain of the ice sheet, H, is divided into a subdomain (or a 
collection of subdomains) ^siA, in which the SIA is sufficiently accurate, and 
in the complement flps = ^\^siA where the FS equations must be solved, 
see Fig. 3. The nodes in the solution vector are reordered according to which 
domain they belong. Then (11) can be rewritten as a 2x 2-block system, one part 
representing the FS subdomains, and one part representing the SIA subdomains 


Au 


a _ 
FS — 


f ^FS 
yAoc 


Aco\ 

Asia] 


/ bfifS \ 

J ■ 


( 12 ) 


In this equation, Aps is the FS system matrix for the FS area, Asia corresponds 
to the SIA area, Aco represents the dependency of the solution in the FS area 
on the values of u in the SIA area, and Aqc represents the dependency of the 
solution in the SIA area on the solution in the FS area. Both Aps and Asia are 
square matrices. In the beginning of each time-step, usia is computed for the 
whole domain O by evaluating the expressions in (9). This computational work 
is negligible in comparison to solving the FS equations. The solution can 

then be approximated by in flsiA, and we only have to solve for 




Figure 3: The domain n is automatically partitioned into subdomains Sds/A 
where the SIA is sufficiently accurate, and other subdomains ilps where the FS 
equations need to be solved such that Q — V,fs^ ^sia- 


disregarding the entire second row in (12). The remaining equation in the top 
row of (12) is then 

+ AcougzT = . (13) 

The solution is an approximation of since is replaced by 

The SIA solution is known in fl and we can move the second 

term to the right hand side in (13) to arrive at an equation for u^g® 

Afsu^^ = b""®® - AcougfA". (14) 

In this way, the SIA solution will act as Dirichlet boundary condi¬ 

tions for the FS solution at the interfaces between Hfs and flsiA- Remember 
that ice is a non-Newtonian fluid with a viscosity dependent on v such that 
Afs = AFs(upg®) and Aqo = Aco(ufs^)- "bhe system (14) has to be solved 
iteratively by a non-linear Picard or Newton iteration, but u^^® has much fewer 
unknowns than u^g in (12). At convergence, the complete ISCAL solution is 
merged together as uJg^AL = (u?scal)^ = ((ufs®)'^> (^siaV)- 

An outline of the algorithm is found in Fig. 4b. In comparison to the FS 
solver in Fig. 4a, error estimation, a SIA solver, and administration of the 
partitioning of 11 are added in the ISCAL solver. 

3.4 Model adaptivity and error estimation 

It is difficult to know a priori exactly where in an ice sheet the SIA is a valid 
approximation, and the answer will change when the ice sheet evolves in time. 
Therefore, we use an automatic error estimation of a linearized problem, which 
determines in which parts of the ice sheet the SIA can be applied, and split the 
domain 11 into ilsiA and Hfs accordingly. The highest update frequency for 
the splitting of the domain is at most at the end of every time step, but for 
most problems this update can be much less frequent. Since there are several 
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Initial guess 
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-i 

Update ice surface 


(a) FS 


Initial guess u^ppg, 
initial division il — Qpg U ^^siA 

I 

- > time-stepping 


Compute 


^ Non-linear 
iteration 


Compute T] = 

Assemble Afs andApo 

I 

Solve 

AFs4f = h^^s.Acon^/;-/ 

I 

11 ^ — 


4 


Error ^ 
estimation ^ 

I 

No 

1 

Update ice surface ^ 

-i 

(b) ISCAL 


Assemble 

■A-F5(uffs) 

+ 

Compute 
error estimation 

Sort nodes into 
fVs and UsiA 


Figure 4: Algorithms for computing ice sheet evolution by solving the FS equa¬ 
tions and the ISCAL equations. 
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possible measures of the accuracy of the SIA, we construct three different error 
estimations; one based on the velocity, one on the residual of (2), and one on 
a functional of the velocity. The parameters for the error estimations, such as 
the error tolerance, e, and the number of timesteps there are between each error 
estimation, to, are user defined. 


3.4.1 Estimating the error in the velocity field 

Determining accurate ice velocities, e.g. for estimating mass loss or for com¬ 
parison with satellite data, is an important glaciological problem. From this 
viewpoint, it is of interest to control the error in the SIA velocity, i.e. in the 
velocity components of 

= UfS ~ (15) 

Since we do not have access to ups, the error is approximated by comparing to 
a reference solution, Uppp, instead, i.e. 


^^SiA = u' 


REF 


11 ^ 


(16) 


The reference solution, uppF, is obtained once a solution u^scal been 
found, by solving the linear system 

= b, (17) 

where A = A(upg.(;;.^j^). For comparison remember that the FS solution satisfies 

A(u9_5)u^ 5 = b. (18) 

In order for the error estimation to be accurate, i.e. for « Au§j^, the 

difference between Upg and should be small. By denoting the error in 

^%CAL by Au/sc'.4l such that 

Au/scAL = UpF — UpgFALJ (19) 


and using (18) we get 

M^jscal + Au/scal)(uj5C’al + Au/scal) 

= (A -h A(up_g^^^ AuisCAl) - M^?SCAL))i^?SCAL + ^^iscal) (20) 
= Au9s + dAAuiscALUps = b, 

where an element dAiji^ of 5A is the first derivative of with respect to u^. 
From (20) we conclude that if ||(9 AAu/sc’AL Upg|| is small then in (17) 

is a good approximation of in (19) and Au^j^ « Au^j^. 

The computation of only requires the assembly and solution of a 

system of linear equations. The extra computational work is acceptable, in 
particular since the error estimation is generally not necessary in every time 
step. The additional computational time attributed to the error estimation is 
measured in Section 4.2. 
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In an ice sheet, the major, most easily observed flow of ice is in the hori¬ 
zontal plane, and we therefore control the error only in the horizontal velocities 
(vx,Vy) = {uv^,Uvy). It is, however, easy to adapt the code into estimating the 
error in all components of the solution. To determine whether SIA is sufficiently 
accurate in a node i in the domain, both the relative and absolute differences 
between the nodal values (vxi,Vyi)siA and (vxi,Vyi)REF are considered at each 
node {i = 1, ... ,N), such that if 

II = \\{Vxi,Vyi)siA — {Vxi,Vyi)REF\\ 

< Taayi{erel\\(Vxi,Vyi)REF\\,iabs) 

then the SIA is sufficiently accurate, and the node is included in ^siA- The 
absolute and relative error tolerances, eabs and Creh are user defined. In the right 
hand side of (21), the absolute error tolerance is chosen when \\{vxi,Vyi)REF\\ 
is very small and a low relative error may not be necessary. Otherwise, the 
relative error tolerance is chosen. 


3.4.2 Estimating the error in the residual 

A computationally cheaper approach is to consider the error in the residual. 
It measures how well the discretized equations are satisfied by the numerical 
solution. The residual of the equations in (2a) and (2b) is a right hand side r 
which is non-zero if there is an error in the solution. 

Hence, the residual in the numerical solution of (11) should be as small as 
possible. 

Again let Auria be the deviation of from the reference solution 
in (17). Then the residual for — Aus/^ can be estimated as 

rsiA = b — = A{u^pp — Ugj^) = —AAus/^. (22) 


Only a matrix-vector multiplication is required to compute tria- Again, if 
A = A, the residual would be exact. The residual corresponding to node i is 
denoted by tria^i- If |lfs'/A,i|| > Cres for a given tolerance Cres, then this node 
is included in V,er- The difference between unreal non-zero only 

in flFS- Since 

=b^--(Ap5ugfi+Acoug,T) 

= {Afru^^ + Acougfi") - {AFRU^fX + Acougfi") (23) 

= AFR{n%^-n2fX), 


the error in the solution and the error in the residual are related as u 


FS _ A-1 


^SIA 


FS 

SIA- 


Ofs 

FS 


3.4.3 Estimating the error in a functional of the solution 

In glaciology, sometimes the velocity field is not of primary interest, but rather 
some functional of the velocity, such as the ice discharge from a drainage basin 
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(see [16]). Let /(u) be a linear functional of the solution with /(O) = 0. On 
the discrete mesh, this functional can be written as f^u, i.e. the scalar product 
of a vector f defining the functional on the mesh and the solution vector u. To 
understand which parts of the ice sheet are important for this functional, the 
adjoint equation or dual problem 


A^w = f (24) 

is solved for the weight vector, w, expressing the impact on / of the solution u 
in each mesh node. We solve the linear system (24) with A = A(u)(g^^j^) as 
in Section 3.4.1. Transposing (24), multiplying with the solution error Au, and 
inserting the residual yields the error in the functional f^u as 

f^Au = w^AAu = w^r. (25) 

Note that in this case we can either compute f^Au directly, or evaulate w^r. 
We choose the latter. The residual is evaluated as in Section 3.4.2, meaning 
that r = isiA and Au = Aug/^. Let N be the total number of nodes in O. If 
Wi is the weight vector corresponding to node i and \'wfrsiA,i\ > e//-^) then 
the error in the functional due to node i is larger than acceptable based on 
equidistribution of the error and the user defined error tolerance e/, and this 
node is included in V,fs- H |wffg/A,i| < e//-^ for all i, then by (25) 

N 

[FA'us/aI = Iw'^fg/^l < \^frsiA,i\ < e/. (26) 

i=l 

This way both the model error in SIA measured by the residual tsia and the 
importance of the node are considered. This approach can be applied to 
any linear functional /(u). It is the discrete version of the goal oriented or 
a posteriori error control in finite elements where the solution to the adjoint 
differential equation provides the weights w [Ilj. In Section 3.4.2, the nodal 
weight Wi is equal to YsiA,i- 

3.5 Computational performance 

The computational work for both ISCAL and FS is dominated by solving the 
Stokes problem. The additional cost of the solution of the free surface problem 
in (1), and other auxiliary processes, can be assumed to be small in compari¬ 
son. The main computational work when solving any finite element problem is 
attributed to the assembly of the system matrix and solution of the associated 
linear system. Often, the solution time dominates the assembly time for large 
problem sizes, but in some cases the assembly time is significant. For instance 
in non-linear problems, the solution time is often reduced by using the solution 
from a previous non-linear iteration as an initial guess for the linear system 
solver, while the assembly time remains large. When solving the FS system in 
(II), the entire system matrix. A, is assembled. In ISCAL, only the smaller 
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matrices Aps and Aqo are assembled, which is less costly. The solution time 
is also lower for ISCAL, since (14) is smaller than the original system (11). 

Let Nps be the number of nodes in ftps, and Nsia = N — Nps be the 
number of nodes in ^siA- Assume that the number of non-linear Picard or 
Newton iterations required to solve (11) and (14) is ^ps and Fiscal, respec¬ 
tively. Furthermore, assume that the work to solve the linear systems once is 
CpN^ and CsNpg for FS and ISCAL respectively and that the work for the 
assembly of the system matrices (A, or Aps) is CaN^ and CaN^q. Both Cp 
and Ca are constants. The choice of linear solver determines Cp and S, and the 
implementation of the assembly determines Ca and A. In ISCAL, there is extra 
work included in the assembly phase because not only is the system matrix, 
Aps, assembled but also Aco, and because of the reorganization of (II) into 
(14). We represent this extra load by adding a constant Co to Ca- Assuming 
that any process beside solving the Stokes system is negligible, the speedup, q, 
of ISCAL can be described by, 

FS load 
ISCAL load 

^_ ^fs{CaN^ + CsN^) _ 

^iscal{{Ca + Co)N^s + CsN^s + ' 

If the system is solved by Gaussian elimination then S' = 3 and if it is solved 
by an optimal multigrid algorithm S may approach 1 [7]. The system matrix is 
sparse, and so A = 1. 

The term {CaN^ + CsN^)/{m^iscAL) represents the error estimation. It 
is small because of the factor {m^jscAL), be. because the error estimation 
is not done in every non-linear iteration, and because the error estimation is 
typically not done in every timestep, i.e. m > 1. If the residual based error 
estimate is used, CpN^ is removed. 

If we assume that Fiscal is almost equal to ^ps, and that Co is not much 
greater than Ca, then we have q > 1, since Nps < N. In Section 4.2, the 
computational gain of using ISCAL is measured in numerical experiments and 
(27) is revisited. 

4 Numerical Experiments 

In this section four experiments are performed. In Experiment I, the ISCAL 
is evaluated for accuracy and efficiency. The ability to dynamically change the 
FS domain Vlps is tested in Experiment 2. In Experiment 3, the three different 
error estimates are compared. ISCAL is tested on real data in Experiment 4, . 
In Experiments 1-3 a conceptual model problem (see Section 4.1) is considered, 
and in Experiment 4 data from the Greenland Ice Sheet is used. The simulations 
are performed on a single core, although a parallel version of ISCAL now exists 
and is used in [24]. 
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4.1 Model problem description 

The geometry of the first three experiments is a 3D circular ice sheet. The 
upper surface position, h{x,y,t), is initialized as a so-called Vialov prohle [17] 

/ („+l)/„\ "/(2"+2) 

h(a;,y,0) = ho fl-(-) j + 100, (28) 

and evolves in time according to (1). Here /iq = 3575.1 m is the maximum 
initial height, L = 750 km is the radius of the ice sheet, and r is the distance of 
a point (x, y) from the center. By comparison, the Greenland Ice Sheet is about 
2400 km long and 1100 km wide, and is about 2000 to 3000 m thick. 

The bedrock is flat {i.e. 6 = 0), and since the problem is isothermal, the rate 
factor is constant, A — 10“^® Pa~^ year“^. The Glen parameter, n, is set to the 
standard value 3, the density, p, is 910kgm“^, and the accumulation/ablation 
function is given by 

=min{0.5,10-^(450 • 10^-y/x2+p^)} (29) 

m/year, so that the the accumulation is positive in the center of the ice sheet 
and decreases radially, becoming negative at r = 450 km, mimicking accumula¬ 
tion/ablation patterns observed in natural ice sheets. Note that this is not the 
accumulation leading to the Vialov profile in [17], but rather the accumulation 
function used in the EISMINT benchmark experiment [21]. At the base of the 
ice, a no-slip condition is applied for Experiment 1 and 3, and a linear slip 
condition for Experiment 2. 

The mesh is constructed by vertically extruding a 2D footprint mesh to 20 
layers (see Fig. 2). The 2D footprint mesh is unstructured and triangular in 
the center of the circle, and radially aligned towards the margins. 

In all experiments, the linear system for the Stokes problem is solved it¬ 
eratively with the generalized conjugate residual method, and the non-linear 
equation is solved with a Picard iteration. An under-relaxation with factor 0.9 
is applied to the Picard iteration in order to stabilize it. 

4.2 Experiment 1 - Efficiency and accuracy 

4.2.1 Setup 

The purpose of ISGAL is to reduce the simulation time with a controlled, small 
reduction in accuracy. In this experiment, four different meshes of varying 
resolution are employed to simulate the model problem with ISCAL, FS, and 
SIA in order to compare simulation time and accuracy. The coarsest mesh has 
17860 nodes (see Fig. 2), and the finest mesh has 257000 nodes. The number 
of nodes in the finest mesh is of the same order of magnitude as the number of 
nodes used for simulating the Greenland Ice Sheet in [35] and [16]. Since the 
number of non-linear iterations varies in time, we run the simulation for 30 time- 
steps of one month each. When measuring the CPU-time, the un-representative 
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Figure 5: Experiment 1 - The velocity magnitude of a circular ice sheet, as seen 
from above, after 30 months. The solution is computed with ISCAL on a mesh 
with 96520 nodes. 


first time-step is excluded. In ISCAL, we estimate the error in the horizontal 
velocities every tenth time step as in (21), with the relative tolerance, Crei, set 
to 5% and the absolute tolerance eats = 1 m/year. 

4.2.2 Results 

The magnitude of the velocity is high at steep margins (1000 m/year or more), 
and low near the dome, (Fig. 5). The same pattern is observed in nature, e.g. 
in the Greenland Ice Sheet [22]. The velocity field in Fig. 5 is computed by 
ISCAL. The FS and SIA solutions are not shown since the difference between 
ISCAL, FS, and SIA velocity is difficult to discern by the naked eye (using a 
reasonable color scale), despite considerable relative errors in the SIA (Fig. 6a). 
As expected from theory {e.g. [5, 20, 2]), the SIA error is particularly high 
near the margins and at the dome (in some places greater than 100 % ) (Fig. 
6a). ISCAL automatically detects these problematic regions through the error 
estimation, and applies the FS where the error is higher than 5 % or 1 m/year 
(Fig. 7a). The estimated error from (16), shown in Fig. 6b, is very close to the 
true error from (15) and hence is a good base for the error control. Note that 
the linkage between the SIA and the FS is very smooth, and not visible in Fig. 
5. The relative error of ISCAL compared to the solution of the FS equations is 
much lower than for the SIA (Fig. 7b). In this case, the FS is not applied at the 
dome even if the relative error is high there. This is so because the horizontal 
velocity is very low at the dome and we allow for an absolute horizontal velocity 
error of 1 m/year, see (21). 

The CPU-time of the ISCAL is considerably lower than for the FS system 
(Fig. 8). The speedup increases with problem size, and on the finest mesh 
ISCAL is nine times faster than the FS. The CPU-time required to compute the 
SIA is negligible compared to both FS and ISCAL, which is why it is so popular. 
Table 1 provides more insight in what operations are the most time-consuming. 
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(a) Relative Error in SIA, Ausia (b) Estimated Error, Ausia 


Figure 6: Experiment 1 - Exact and estimated relative error in SIA in percent, 
after 30 months, viewed from above. The relative error tolerance, ereh in SIA 
is 5 % and the absolute error tolerance, Cats is I m/year. The color scale is set 
so that everything less than 5 % is blue. 




% 

^00 

lio 

I 

1 


(a) ^ps (red) and ^^siA (blue) (b) Relative Error in ISCAL, AujgcAL 

Figure 7: Experiment 1 - The distribution of flps (red) and flsiA (blue) and 
the error in ISCAL compared to FS after 30 time-steps. 
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Figure 8: Experiment 1 - CPU-time versus number of nodes for the SIA (blue 
dashed line), ISCAL (black solid line) and the full Stokes (red dashed line) 
problem. The labels indicate speedup, q, of ISCAL relative to FS. 


for both FS and ISCAL, on all four meshes. For FS, the assembly and solution 
of the Stokes system completely dominates other processes, such as calculation 
of the SIA or solving for the free surface position, as assumed in (27). For 
ISCAL, also the error estimation time is significant, although still much smaller 
than the time for the assembly and solution. The relative number of FS nodes 
decreases with problem size, most likely because the numerical errors in both 
SIA and FS are reduced with increasing resolution. Because of this, the relative 
time spent on the error estimation increases with the problem size, whereas 
the assembly and solution time decreases, while the error estimation time is 
independent from the size of ^.ps- The assembly time dominates the solution 
time, even though the importance of the solution time increases slightly with 
problem size as expected if S' > 1. The average number of non-linear iterations 
per time step increases slightly with problem size, and is the same for ISCAL and 
FS, because the areas that converge slowly are included in ^fs- By revisiting 
(27), using the data in Table 1, we can better understand how the speedup, q, 
is related to N/Nps- If we use that 1) Ca > Cs in (27), 2) S is not very much 
larger than unity, 3) the time spent on error estimation is small, and 4) that 
^FS ~ Fiscal, we can roughly approximate (27) by 

q = ^ ^ A ——f small remainder. (30) 

Ca + Co Nfs 

Fitting q{N/NFs) to a linear polynomial gives Ca/{Ca + Co) ~ 0.53, such that 
the speedup is less than N/Nps (Fig. 9). In other words, the computational 
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Table 1: Experiment 1 - Percentage of full Stokes nodes, percentage of CPU¬ 
time spent on matrix assembly, linear system solution, and error calculation, 
and average number of non-linear iterations per time step for the four different 
meshes. 


# nodes, N Model 

FS nodes 
(%) 

Assembly 

(%) 

Solve 

(%) 

Error Calc. 
(%) 

# iter. 

FS 

100.0 

87.0 

12.3 

- 

10.4 

17860 ISCAL 

26.8 

84.8 

10.7 

2.1 

10.2 

FS 

100.0 

86.6 

12.7 

- 

11.6 

96520 ISCAL 

12.2 

82.3 

9.6 

3.8 

11.9 

FS 

100.0 

85.7 

13.7 

- 

12.7 

175960 ISCAL 

9.1 

82.0 

00 

bo 

4.4 

12.7 

FS 

100.0 

84.8 

14.6 

- 

12.9 

257000 ISCAL 

6.8 

78.3 

11.2 

5.1 

12.9 


work associated with rearranging (11) to (14), and the assembly of Aco, is of the 
same order of magnitude as the assembly of A.fs in the present implementation. 
This reduces the speedup, q. In a general context, often Ca < Cs, which would 
make the speedup larger. The speedup would also be larger if a direct solver is 
used, since then S' « 3. As indicated in Fig. 9, increased speedup with problem 
size can be attributed to the fact that the proportion of FS nodes decreases with 
increasing problem size. There are, however, additional explanations, neglected 
in (30), e.g. that the overhead from sorting the nodes has a larger relative 
impact for small problem sizes, and that S > 1. 

4.3 Experiment 2 - Ice stream tracking 

4.3.1 Setup 

In this experiment, we test how well ISCAL adapts to rapid changes in dynamics, 
and how it treats fast flowing ice streams. Ice streams are very important for 
the mass budget of the Greenland Ice Sheet, and even more so for the Antarctic 
Ice Sheet [33, 23]. As the SI A assumes that shearing motion dominates sliding, 
we expect ISCAL to apply the Stokes equations in ice streams. Ice streams may 
change position during long time spans, [10, 41, 9] and ^Ips should adapt to 
this change. In order to test this, the basal boundary condition is changed from 
no slip to the linear sliding law in (8) with 


P{x, y, t) = max 




(31) 


where 0{x, y) is the polar angle of every point (x, y). This sliding law allows for 
an ice stream to form, starting at the center and flowing radially outwards to 
the margin, moving counter clockwise as the time t increases. The area outside 
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Figure 9: Experiment 1 - Speedup q versus percentage of FS nodes (stars). 
The dashed lines indicates the polynomial fit g = 0.53{N/Nps) + 1.47. In our 
experiments, the larger the problem size, the smaller the percentage of FS nodes. 


the ice stream has a very high friction (/3 = 10“*^'^) effectively imposing a no 
slip condition. ISCAL is applied using the error estimation based on horizontal 
velocity (Section 3.4.1) with a relative tolerance, e^eh of 5 % and eabs = 1 
m/year, on a mesh with 96520 nodes. The error is estimated in every time-step. 
If the error estimation is performed less frequently, there will be a delay in the 
update of flps and ^IsiA- 

4.3.2 Results 

The modified sliding law creates an ice stream in which the velocity is signifi¬ 
cantly higher than in the surrounding ice (Fig. 10). ISCAL automatically tracks 
this ice stream as it moves, applying the Stokes equations in the fast flowing 
region and at the margins as expected (Fig. 11). 

4.4 Experiment 3 - Comparing error estimations 

4.4.1 Setup 

The three different methods of error estimation presented in Section 3.4 will lead 
to different distributions of the FS nodes, i.e. of ilps and ilsiA- To compare 
the different distributions we run three simulations using each one of the three 
estimates. The set-up is the same as in Experiment 1, and the simulations 
run for 12 time-steps (months) using the mesh with 96520 nodes. We keep the 
number of nodes in ftps constant at 20 % of the total {i.e. 19304 nodes), such 
that the 19304 nodes with the highest error, computed with respective error 
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Figure 10: Experiment 2 - Velocity magnitude with ISCAL after 11 months, 
seen from above. The boundary between Vlps and ^siA is not visible. 



(a) 1 month (b) 6 months (c) 11 months 

Figure 11: Experiment 2 - Distribution of Hfs (red) and ^siA (blue) after 1, 
6 and 11 months. The ice stream position switches counter clockwise. 
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(a) Horizontal Velocity 


(b) Flux 



(c) Residual 


Figure 12: Experiment 3 - Distribution of ftps (red) and flsiA (blue) for the 
three different error estimations after 12 months. A slice is cut out to show how 
the FS nodes are distributed within the ice, and the vertical axis is scaled by 
a factor 100. The vertical black lines in (b) indicate the position of the surface 
F : = 600 km. 

measures, are included in Ofs- When estimating the error in the velocity field, 
the nodes are sorted based on the relative error in the velocity, but a node that 
has an absolute error that is lower than an absolute tolerance of 0.1 m/year 
will not be included in Hfs- This prevents areas with almost zero velocity {e.g. 
near the base) from inclusion in Hfs- The absolute tolerance is was lowered 
compared to the previous experiments, in order to avoid repetition. For the 
functional based error estimation, (Section 3.4.3), the functional is chosen as 
the ice flux across the surface F{x, y) : = 600 km 



(32) 


where ny is the normal pointing outwards from the surface F. In a real ice 
sheet, a large part of the mass loss is through ocean terminating ice streams. A 
functional describing the ice flux across a vertical surface can be used to control 
the error in the discharge through such streams. The simulation where the error 
in the residual is estimated is straightforward. 

4.4.2 Results 

As expected, the FS area, Hfs, depends on the method of error estimation (Fig. 
12). When estimating the error in the horizontal velocity, the distribution is 
similar to the previous experiments, i.e. FS is applied at the margins (Fig. 12a), 
but as the absolute tolerance was lowered from 1 m/year to 0.1 m/year also the 
dome is now included in SIfs- When estimating the error in the ice flux, the FS 
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nodes are concentrated around the surface F{x,y) (Fig. 12b). As there is no 
sliding at the base, the main motion is by vertical shearing, such that the flux 
through the surface is higher near the surface than close to the base. Therefore 
there will be more FS nodes near the surface than near the base, as observed 
in Fig. 12b. The simulation estimating the error in the residual results in FS 
nodes at the dome, near the margins, and close to the ice surface (Fig. 12c). 
Since the solution is numerically sensitive near the surface, the distribution of 
nodes is somewhat irregular there. The same effect is found in Fig. 12a. The 
magnitude of the error is higher at the margins than at the dome and at the 
surface. Thus if the number of FS nodes is decreased, ^fs is similar for the 
residual based and the velocity based error estimates. The computational work 
to estimate the error in the residual is lower than in the other error estimates, 
because no linear system is solved. On the other hand it is perhaps not as 
intuitive as estimating the error in the velocity or in the flux. 

4.5 Experiment 4 - Application to the Greenland Ice Sheet 

4.5.1 Setup 

In order to test ISCAL on realistic data of interest to the glaciological commu¬ 
nity, we apply the ISCAL to the Greenland Ice Sheet. Note that we are not 
aiming to obtain results that give insight into the ice sheet dynamics of Green¬ 
land nor to answer any glaciological questions. We would like to emphasize that 
we simply want to use the Greenland data in order to test the method on a real- 
world geometry, rather than to present a complete simulation of the Greenland 
Ice Sheet. The application of the method to problems of glaciological interest 
will be reported in future studies. For simplicity, we consider isothermal con¬ 
ditions (with A as in previous experiments), do not apply any time-integration 
and thus need no climatic forcing (z.e. = 0). The data required are bedrock 

topography, &, under the ice sheet, ice surface topography, /i, and a sliding co¬ 
efficient P describing the friction between the ice and the bedrock. The sliding 
coefficient is incorporated in the linear sliding law of (8). 

Both the bedrock and the surface topography (Fig. 13) are from measure¬ 
ment data in [4], while the sliding coefficient /3 in (8) was computed by inverse 
modeling in [16]. Gontrary to the model problem, the bedrock is far from flat, 
with the East Greenland mountain range, and large areas below sea level in 
the interior (Fig. 13b). The sliding coefficient is high {i.e. high friction) in 
the south, and low in many outlet glaciers, e.g. in the large ice stream in the 
northeast, NEGIS (Northeast Greenland Ice Stream) (Eig. 13c). 

As the SIA is sensitive to high frequency spatial variations [19, 2], the data 
are smoothed to a resolution of 50 km (which is what is shown in Eig. 13). 
We generate a triangular mesh with an edge length of 10-20 km. To simplify 
the mesh construction and avoid small skewed elements at bays and capes, 
we smooth the margin slightly. The error estimation is based on the relative 
horizontal velocities, using a relative error tolerance as in Experiment 1 and 2, 
i.e. Crei is 5 % and an absolute error tolerance. Cats, of 1 m/year. 
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(a) Surface topography, h (b) Basal topography, b (c) Basal sliding coeffi¬ 

cient, /3 

Figure 13: Experiment 4 - Input data for the simulation of the Greenland ice 
sheet. The geometry is from meausured data from [4], and the sliding coefficient 
is computed in [16]. The data is smoothened to a resolution of 50 km. 


4.5.2 Results 

The computed ISCAL velocity field is shown in Fig. 14a, and compared with 
observational surface (horizontal) velocity data from [22] (Fig. 14b). As with 
the model problem, the velocity is high at the margins and in places with low 
friction, and low at the domes. The general pattern of the ISCAL velocity 
matches observations but the flow velocity is in general too high. The reason 
is because we did not allow the errors in the prescribed initial ice geometry 
to vanish, and since the data was smoothened and the mesh is fairly uniform. 
When simulating a real ice sheet, the velocity field and free surface are is usually 
allowed to relax, great care is taken into the initialization of the simulation. 
A realistic, transient, fully resolved, simulation of the Greenland Ice Sheet is 
however beyond the scope of this paper. The resolution is too low to capture 
all the small outlet glaciers at the margins, but large features such as the ice 
stream NEGIS in the northeast, the Petermann glacier in the northwest, and 
Jakobshavn Isbrae on the west coast are clearly visible. The error in SIA in Fig. 
15a is high at the margins and domes (here ridges), cf. the model problem in 
Section 4.2. Also in the high sliding area under the NEGIS in the northeast, the 
error is high. For this diagnostic simulation, the estimated SIA error and the 
true SIA error will be equal, as we always start a simulation with flps = ^- The 
distribution of and Hs/a in Fig. 15b follows from the error. As expected 
when using real-world data, the pattern is more irregular than for the model 
problem, but the general behavior is the same. The error tolerances are very 
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(a) Velocity, ISCAL (b) Observed velocity 


Figure 14: Experiment 4 ISCAL velocity and observed velocity (re-plotted from 
[22]) and displayed on a common colorscale. White patches in the observed 
velocity are due to missing data. Three ice streams/outlet glaciers discussed in 
the text are named. 
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(a) Relative Error in SIA, (b) flps (red) and QsiA 
Ausia (blue) 

Figure 15: Experiment 4 - Estimated relative error of horizontal velocities in 
percent, and distribution of flps and flsiA- The relative error tolerance, erei, 
in SIA is 5 % and the absolute error tolerance, eats is 1 m/year. 
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low, and considering the uncertainty in data they can be allowed to be higher. 
If a relaxation of the velocity held and surface would have been applied, we 
expect the error and the extent of ^Ips to decrease. 


5 Conclusions 

ISCAL drastically reduces simulation time while keeping errors caused by the 
approximation of the Stokes equations at acceptable and controllable levels, by 
solving the full Stokes (FS) equations only where needed and using the SIA 
elsewhere. For a problem size of the same order of magnitude as for realistic 
ice sheet simulations, ISCAL is nine times as fast as FS, while the error in the 
horizontal velocity is no greater than 5 % or I m/year. The speedup increases 
with problem size, mainly because the relative number of FS nodes is decreasing 
with increasing resolution. ISCAL automatically detects rapid changes in ice 
dynamics and applies the FS equations in regions where the SIA is not accurate 
enough, such as at margins, in ice streams, and at domes. The estimated error 
is close to the true error and also agrees with the theory of the SIA. Three 
different ways of determining the partial domain where to apply the SIA are 
implemented; estimating the error in the horizontal velocity, estimating the 
error through a functional of the velocity (ice flux), and estimating the error 
utilizing the residual of the Stokes equations. Each error estimate results in 
different distributions of the areas where the FS equations are solved, and are 
suitable for different applications. We further show that ISCAL is applicable 
to real problems, such as the Greenland Ice Sheet. The error and distribution 
of FS areas are more fragmented when using real, less smooth data, but the 
general pattern is the same and the SIA equations are sufficiently accurate in 
large parts of the ice sheet. 

When using the ISCAL method, special care should be given to the meshing 
procedure. As the SIA is sensitive to high frequency spatial variation in data 
(e.g. roughness in underlying bed), the data might have to be smoothened in 
order to avoid unnecessarily high SIA errors. Also, if the mesh is extra refined 
in areas where the SIA will likely not be valid, such as in ice streams, the 
speedup of the ISCAL will be lower. The ISCAL method is developed with 
paleosimulations in mind, i.e. large scale ice sheet simulations over long time 
spans. As the uncertainties in paleosimulations are anyway very large because 
of lack of accurate data, relatively coarse grids and smoothened data are not a 
problem. 

Due to the computational costs of using the exact FS equations, the SIA 
and other approximations have long been the only option for paleosimulations, 
despite high errors in crucial areas. Considering how complex and variable the 
dynamics of a natural ice sheet is, it is obvious that any approximation will 
result in significant, over time accumulated, errors in some parts of the domain. 
We believe it is important not to apply approximations without controlling the 
modeling error, and ISCAL is a powerful tool for this. 
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